##------------------------------------------##
## Lauren Peritz
## ISQ 2019 - WTO Compliance
## Replication Files - Simulations for Hypothetical Trade figure 1                
## Update 2019.10.30
##------------------------------------------##

rm(list=ls(all=TRUE))


set.seed(12345)

#### Make some trended graph that looks like trade 
time <- seq(1,20,1)
Y.observed <-  0.6 + .1*runif(20)-0.12*log(time)
Y.counterfactual1 <-  c(Y.observed[1:10]-.01+0.025*runif(5), Y.observed[11:20]+0.5*seq(0.01,0.1,.01)+0.04*runif(10))
Y.counterfactual2 <-  c(Y.observed[1:10]-.01+0.025*runif(5), Y.observed[11:20]-1.5*seq(0.01,0.1,.01)-0.08*runif(10))
Y.counterfactual3 <-  c(Y.observed[1:10]-.01+0.025*runif(5), Y.observed[11:20]-.02+0.035*runif(5))

#par(mfrow=c(1,2))
par(mar = c(8,4,4,2))

# Set background color
library(scales)
tr.gray <- alpha("gray", alpha=.3)
tr.white <- alpha("white", alpha=.2)

#### Hypothetical Compliance Graph
# Step 1
plot(time, Y.observed, col="white",
	xlim = c(1.5, 19.5),
	ylim = c(.05,.7), "l", lwd=2,
	ylab = "", xlab = "", xaxt="n", yaxt="n",
	main = "(a) Compliance Trade Pattern")

# Step 2
mtext(side=2, line=1, "Trade", cex=1.5)
mtext(side=1, line=7, "Time", cex=1.5)
rect(xleft=6, xright=9.8,ybottom=0, ytop=.75, col=tr.gray, border = tr.gray)
abline(v=9.8, lwd=2, col="grey35")

# Step 3
lines(time, Y.observed, lty=1, lwd=2)
legend("topright", lwd = c(2), lty = c(1), c("Treated"), bty="n")


# Step 4
lines(time, Y.counterfactual2, lty=2, lwd=2)
legend("topright", lwd = c(2,2), lty = c(1,2), c("", "Control"), bty="n")
polygon(c(time[10:20], rev(time[10:20])), c(Y.observed[10:20], rev(Y.counterfactual2[10:20])),
	density=15, angle=90, border=NA, col="grey")
polygon(c(time[10:20], rev(time[10:20])), c(Y.observed[10:20], rev(Y.counterfactual2[10:20])),
	col=tr.white, border=NA)





#### Hypothetical Noncompliance Graph
plot(time, Y.observed, 
	xlim = c(1.5, 19.5),
	ylim = c(.05,.7), "l", lwd=2,
	ylab = "", xlab = "", xaxt="n", yaxt="n",
	main = "(b) Noncompliance Trade Pattern")
lines(time, Y.counterfactual1, lty=2, lwd=2)
mtext(side=2, line=1, "Trade", cex=1.5)
mtext(side=1, line=7, "Time", cex=1.5)
legend("topright", lwd = c(2,2), lty = c(1,2), c("Treated", "Control"), bty="n")

## Indicate dispute period and treatment time
rect(xleft=6, xright=9.8,ybottom=0, ytop=.75, col=tr.gray, border = tr.gray)
abline(v=9.8, lwd=2, col="grey35")

## polygon fill
polygon(c(time[10:20], rev(time[10:20])), 
	c(Y.counterfactual1[10:20], rev(Y.observed[10:20])),
	density=15, angle=90, border=NA, col="grey")
polygon(c(time[10:20], rev(time[10:20])), 
	c(Y.counterfactual1[10:20], rev(Y.observed[10:20])),
	col=tr.white, border=NA)



























##################################### old version #####################







#### Make some trended graph that looks like trade 
time <- seq(1,20,1)
Y.observed <-  0.6 + .1*runif(20)-0.12*log(time)
Y.counterfactual1 <-  c(Y.observed[1:10], Y.observed[11:20]+1.5*seq(0.01,0.1,.01)+0.08*runif(10))
Y.counterfactual2 <-  c(Y.observed[1:10], Y.observed[11:20]-1.5*seq(0.01,0.1,.01)-0.08*runif(10))

par(mfrow=c(1,2))
par(mar = c(8,4,4,2))

# Set background color
library(scales)
tr.gray <- alpha("gray", alpha=.3)
tr.white <- alpha("white", alpha=.7)

#### Hypothetical Compliance Graph
plot(time, Y.observed, 
	xlim = c(1.5, 19.5),
	ylim = c(.05,.7), "l", lwd=2,
	ylab = "", xlab = "", xaxt="n", yaxt="n",
	main = "(a) Compliance Trade Pattern")
lines(time, Y.counterfactual2, lty=2, lwd=2)
mtext(side=2, line=1, "Trade", cex=1.5)
mtext(side=1, line=7, "Time", cex=1.5)
legend("topright", lwd = c(2,2), lty = c(1,2), c("Treated", "Control"), bty="n")

## Indicate dispute period and treatment time
rect(xleft=6, xright=9.8,ybottom=0, ytop=.75, col=tr.gray, border = tr.gray)
abline(v=9.8, lwd=2, col="grey35")

## polygon fill
polygon(c(time[10:20], rev(time[10:20])), c(Y.observed[10:20], rev(Y.counterfactual2[10:20])),
	density=15, angle=90, border=NA)
polygon(c(time[10:20], rev(time[10:20])), c(Y.observed[10:20], rev(Y.counterfactual2[10:20])),
	col=tr.white, border=NA)

#### Hypothetical Noncompliance Graph
plot(time, Y.observed, 
	xlim = c(1.5, 19.5),
	ylim = c(.05,.7), "l", lwd=2,
	ylab = "", xlab = "", xaxt="n", yaxt="n",
	main = "(b) Noncompliance Trade Pattern")
lines(time, Y.counterfactual1, lty=2, lwd=2)
mtext(side=2, line=1, "Trade", cex=1.5)
mtext(side=1, line=7, "Time", cex=1.5)
legend("topright", lwd = c(2,2), lty = c(1,2), c("Treated", "Control"), bty="n")

## Indicate dispute period and treatment time
rect(xleft=6, xright=9.8,ybottom=0, ytop=.75, col=tr.gray, border = tr.gray)
abline(v=9.8, lwd=2, col="grey35")

## polygon fill
polygon(c(time[10:20], rev(time[10:20])), 
	c(Y.counterfactual1[10:20], rev(Y.observed[10:20])),
	density=15, angle=90, border=NA)
polygon(c(time[10:20], rev(time[10:20])), 
	c(Y.counterfactual1[10:20], rev(Y.observed[10:20])),
	col=tr.white, border=NA)

